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1.  Introduction 


The  electrical  conductivity  of  materials  is  of  fundamental  interest  in  a  number  of 
research  areas,  many  of  which  are  relevant  to  current  US  Army  Research 
Laboratory  (ARL)  programs.  In  the  Director’s  Strategic  Initiative  (DSI)  “Transport 
in  Complex  Crystalline  Materials  Based  on  van  der  Waals  Heterostructures”, 
heterogeneous  vdW  solids  are  being  fabricated  experimentally  as  possible  next- 
generation  materials  for  electronic  devices.1  Knowledge  of  the  vertical  charge 
transport  occurring  between  layers  is  critical  as  it  may  lead  to  design  and  fabrication 
of  devices  that  are  more  tunable. 

In  another  ARL  application,  the  electrical  conductivity  of  the  boron  carbide  armor 
ceramic,  one  of  the  focus  materials  for  the  Materials  in  Extreme  Dynamic 
Environments  Collaborative  Research  Alliance,  has  been  shown  to  be  a  function  of 
pressure  and  carbon  content.2  Similarly,  the  electrical  conductivity  of  explosive 
product  gases  has  been  shown  to  depend  on  carbon  content,3  and  electrical 
conductivity  models  have  become  a  requirement  for  input  into  continuum-level 
simulations  being  executed  by  researchers  in  the  Multi-Threat  Armor  Branch. 
Although  the  increased  electrical  conductivity  observed  in  detonation  products  of 
condensed  explosives  has  received  considerable  attention  experimentally,4-5 
conductivity  values  corresponding  to  regular  intervals  of  temperature  and  pressure 
are  required  for  accurate  interpolation  in  continuum  simulations. 

In  the  absence  of  experimental  data,  particularly  in  the  case  of  notional  materials  that 
have  not  yet  been  synthesized  (such  as  the  vdW  heterostructures  being  explored  in 
the  aforementioned  DSI),  the  electrical  conductivity  of  materials  can  be  obtained 
using  first-principles  quantum  mechanical  techniques.  One  approach  for  computing 
the  conductivity  from  first  principles  is  based  on  Boltzmann  transport6  theory  where 
the  electrical  conductivity  (a)  is  proportional  to  the  band  velocities  v  and  lifetime  t: 

<*ap  ~  *vp(k)  ,  (1) 

where  a.fi  denote  Cartesian  direction.  The  band  velocities  in  Eq.  1  are  obtained  by 
differentiating  the  band  energies  e  with  respect  to  k  points  in  the  Brillouin  zone — 

ds 

for  example,  v  ~  — .  In  practice,  the  differentiation  is  done  numerically  using  a 

very  dense  grid  of  k-points  as  implemented  in  codes  such  as  BoltzTrap.7  The 
lifetime  r  is  a  parameter  and  can  be  approximated  by  fitting  to  experimental  data, 
which  is  mildly  perplexing  for  materials  where  experimental  data  do  not  exist. 
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Another  approach  for  computing  conductivity  is  based  on  the  Kubo-Greenwood 
(KG)  formalism.8  In  this  approach,  the  frequency  (cu)-dependent  electrical 
conductivity  is  computed  as  a  weighted  sum  over  k-points: 

o-(w)  =  Efc  o'*  (co)  *  W (fc)  ,  (2) 


where  W(k)  is  the  weighting  at  integration  point  k.  ak  (m)  given  by 


os  (of) 


27TB  ih‘L 
3m2(x)£l 


-  F(£J,k)]  x  \(xi’j,^a\%,k)\2S(ej,k  -  ei,k  ~  h<o)  (3) 


where  m  is  the  electron  mass,  Q  is  the  unit  cell  volume,  F  is  the  occupation  number 
for  bands  i  and  j,  e  is  the  band  energy,  and  the  index  a  denotes  Cartesian  direction 
for  the  momentum  operator.  The  KG  approach  has  the  advantage  that  it  does  not 
require  an  estimate  of  the  lifetime  or  numerical  differentiation  of  band  energies, 
and  all  requisite  quantities  can  be  obtained  using  quantum  mechanical  methods 
without  fitting  to  experiment. 

In  this  report,  a  Fortran  90  implementation  of  the  KG  formula,  for  use  with  the 
Quantum  Espresso9  software  package,  is  presented.  The  program,  called  ECON- 
KG,  is  currently  available  at  the  ARL  Defense  Supercomputing  Resource  Center 
but  can  be  easily  ported  to  other  Department  of  Defense  computing  facilities. 
Following  details  of  the  implementation  and  instructions  for  execution,  the 
frequency-dependent  conductivity  of  liquid  sodium  is  presented  and  compared  to 
experimental  results. 


2.  Computational  Methods 


2.1  Implementation 

The  band  energies,  occupation  numbers,  and  momentum  matrix  elements  required 
in  Eq.  3  are  evaluated  using  the  Quantum  Espresso  solid-state  density  functional 
theory  software  package.  However,  the  public  release  version  of  Quantum  Espresso 
only  provides  momentum  elements  where  bands  i  and  j  of  Eq.  3  lie  in  the  valence 
and  conduction  bands,  respectively.  Unfortunately,  this  restriction  is  only 
applicable  for  calculations  with  integer  occupation  numbers  and  for  calculations 
that  use  smearing  (Fermi-Dirac,  Gaussian,  etc.)  where  there  are  non-integer 
occupations  for  bands  near  the  Fermi  level;  matrix  elements  where  indices  i  and  j 
both  lie  in  the  valence  band  also  contribute  to  the  conductivity.  Therefore,  the 
Quantum  Espresso  source  code  was  modified  so  that  momentum  matrix  elements 
between  all  bands  are  now  computed  and  written  to  disk  along  with  the  occupation 
numbers,  band  energies,  and  k-weights  required  to  evaluate  the  conductivity. 
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In  ECON-KG,  the  delta  function  appearing  in  Eq.  3  is  broadened  using  a  Gaussian 
function: 


x  1  v -/,«■  j 

8{£j,k  ~  £i,k  A2  >  (4) 

where  the  broadening  factor  A  is  an  input  parameter  set  by  the  user.  This  parameter 
can  have  significant  effects  on  the  computed  conductivity  as  has  been  extensively 
discussed  by  Knyazev  and  Levashov,10  who  tested  the  sensitivity  of  results  for 
aluminum  using  broadening  factors  ranging  from  0.02  to  0.2  eV.  Pozzo  et  al.11 
suggest  a  factor  that  is  equal  to  the  average  spacing  between  the  computed 
eigenvalues  and  used  values  ranging  from  0.0012  to  0.045  eV  for  their  work  on 
sodium.  In  both  cases,  the  convergence  of  the  conductivity  with  respect  to  the 
broadening  factor  was  improved  by  increasing  the  system  size,  and  all  of  these 
variables  have  to  be  considered  when  converging  values  of  the  conductivity. 

2.2  Execution 

In  order  to  compute  the  conductivity  using  ECON-KG,  the  following  3  steps  are 
necessary: 

1)  Perform  a  single-point  self-consistent  field  (SCF)  calculation  to  get  a  set  of 
valence  and  conduction  bands  for  the  system  using  the  “pw.x”  executable 
from  the  Quantum  Espresso  suite. 

2)  Using  the  converged  bands,  compute  the  momentum  matrix  elements  using 
the  locally  modified  (as  described  previously)  “bands. x”  executable  from 
the  Quantum  Espresso  suite. 

3)  Run  the  ECON-KG  executable. 

Sample  Quantum  Espresso  input  files  (for  steps  1  and  2)  using  liquid  sodium  as  an 
example  are  given  in  the  Appendix. 

Given  a  set  of  bands  and  momentum  matrix  elements,  the  ECON-KG  executable  is 
then  run  to  compute  the  conductivity.  The  ECON-KG  input  file  (“kg.input”)  is 
simply  a  list  of  user-defined  variables  formatted  as  follows: 

13117.18  !  Volume  in  Bohr 

0.05  !  Delta  function  factor  (eV) 

486  !  Number  of  Electrons  in  System 

1.73  !  Fermi  Energy  (eV) 

5.0  !  Fermi  Window  (eV) 

0.001  2.0  0.002  !  Frequency  Range 
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Specifically,  line  1  is  the  unit  cell  volume  in  atomic  units,  line  2  is  the  Gaussian 
broadening  factor  for  Eq.  4  in  eV,  line  3  is  the  number  of  electrons  for  the  system, 
and  line  4  is  the  Fermi  energy  (in  electronvolts)  obtained  from  the  SCF  calculation. 
Fine  5,  the  “Fermi  Window”,  sets  the  range  of  band  energies  (measured  from  the 
Fermi  level)  that  are  included  in  the  calculation.  In  the  input  above,  all  bands  within 
5.0  eV  of  the  Fermi  level  are  included.  As  shown  in  Eq.  3,  the  conductivity  is 
frequency  dependent,  and  ECON-KG  will  provide  conductivity  values  for  all 
frequencies  within  a  range  specified  by  the  user.  The  frequency  range  is  defined  on 
line  6,  in  units  of  electronvolts,  where  the  first  and  second  numbers  are  the  lower 
and  upper  frequency  limits,  respectively,  and  the  third  number  is  the  spacing 
between  points  within  those  bounds.  In  the  example  above,  ECON-KG  will 
compute  the  conductivity  for  values  from  0.001  to  2.0  eV  with  a  spacing  of  0.002 
eV  between  successive  points.  The  conductivity  corresponding  to  zero  frequency 
(“DC”  conductivity)  cannot  be  directly  computed  using  the  KG  expression  since  co 
=  0  introduces  a  singularity.  Therefore,  the  DC  conductivity  must  be  obtained  by 
extrapolation  of  the  low-frequency  values  to  zero. 

2.3  Output 

The  ECON-KG  output  consists  of  a  single  table  where  each  row  contains 
conductivity  values  for  each  frequency  requested  by  the  user.  Each  row  of  the  6- 
column  table  contains  the  following  data: 

Frequency(eV)  Conductivity(S/m)  CTXx(S/m)  CTxx(S/m)  aZz(S/m)  Sum_Rule 

where  column  1  is  the  frequency,  column  2  is  the  trace  of  the  conductivity  tensor 
(Siemens/meter),  and  columns  3,  4,  and  5  are  the  diagonal  components  of  the 
conductivity  tensor  that  contribute  to  the  trace.  The  last  column,  “Sum_Rule”, 
reports  the  value  of  the  following  integral: 

(5) 

where  A  is  the  number  of  electrons  in  the  system.  This  integral,  approximated  using 
the  trapezoid  rule  in  ECON-KG,  should  equal  1  (though  in  practice  it  is  generally 
<1)  and  serves  as  a  check  on  the  quality  of  the  calculation. 

3.  Example 

As  an  example  application,  the  conductivity  of  liquid  sodium  has  been  computed 
using  the  sample  inputs  given  in  the  Appendix.  The  input  structure  was  extracted 
from  a  molecular  dynamics  simulation  of  liquid  sodium,  at  a  temperature  of  950  K, 
using  the  Perdew-Burke-Emzerhof12  density  functional  with  Fermi-Dirac 
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smearing  and  an  electronic  temperature  of  0.082  eV.  A  plot  of  the  optical 
conductivity,  as  a  function  of  k-point  sampling  for  this  structure,  is  shown  in  Fig. 
1.  The  plot  contains  curves  for  multiple  k-point  grids  ranging  from  the  Gamma 
point  up  to  8  x  8  x  8  Monkhorst-Pack13  meshes.  The  experimental  zero  frequency 
value  (2.26  x  106  S/m)  shown  in  the  figure  (red  star)  was  obtained  by  extrapolation 
of  measurements  reported  by  Freedman  and  Robertson14  to  950  K.  As  shown, 
convergence  with  respect  to  k  is  slow,  particularly  in  the  low-frequency  regime. 
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Fig.  1  Optical  conductivity  of  liquid  sodium  using  multiple  k-point  grids.  Experimental 
value  (red  star)  taken  from  Freedman  and  Robertson14. 

The  predicted  DC  conductivities,  obtained  by  linear  extrapolation,  are  given  in 
Table  1,  and  a  plot  of  these  values  versus  k-point  grid  size  is  given  in  Fig.  2.  The 
extrapolations  were  performed  with  the  SciPy  Python  library  using  conductivity 
values  from  0.007  to  0.15  eV  for  each  curve. 

Table  1  DC  conductivity  (106  S/m)  for  each  k-point  grid 


Gamma 

2x2x2 

3x3x3 

4x4x4 

5x5x5 

— • — 

6x6x6 

• 

7x7x7 

# 

8x8x8 

k-Mesh 

a(co=0) 

Gamma 

2.43 

2x2x2 

3.78 

3x3x3 

2.12 

4x4x4 

2.46 

5x5x5 

2.72 

6x6x6 

2.52 

7x7x7 

2.45 

8x8x8 

2.37 

Experiment 

2.26 
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Fig.  2  DC  conductivity  of  liquid  sodium  as  a  function  of  k-point  grid  size 

As  shown  in  Table  1,  the  results  are  in  fair  agreement  with  the  experiment. 
However,  the  error  as  a  function  of  k  is  oscillatory,  and  with  a  mesh  size  of 
6  x  6  x  6,  the  results  are  only  starting  to  show  convergence  (see  Fig.  2). 

4.  Conclusion 


The  example  results  presented  in  this  work  were  obtained  using  a  single  structure 
extracted  from  a  molecular  dynamics  trajectory.  In  practice,  one  should  use 
multiple  structures,  with  appropriate  averaging,  to  determine  the  conductivity 
values.  However,  for  brevity’s  sake,  only  a  single  structure  was  used  for  the 
example  application  presented  in  the  report.  The  results  in  this  work  are  in  fairly 
good  agreement  with  the  experiment.  However,  in  the  author’s  experience,  the  level 
of  agreement  can  shift  dramatically  (factor  of  10  or  more)  depending  on  density 
functional,  system  size,  and  input  parameters  to  ECON-KG.  The  results  can  be  very 
sensitive  to  the  number  of  bands  and  delta  function  width,  and  convergence  with 
respect  to  these  quantities  can  be  challenging. 

Finally,  I  should  acknowledge  Dr  Lazaro  Calderin  of  the  University  of  Florida  who 
I  discovered,  during  development  of  ECON-KG,  was  also  developing  an 
implementation  of  the  KG  formalism  in  conjunction  with  Quantum  Espresso.  I 
thank  Dr  Calderin  for  many  fruitful  discussions  and  for  the  opportunity  to  be  a  beta 
tester  of  his  program15,  which  was  used  to  debug  and  validate  the  results  obtained 
with  ECON-KG. 
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Appendix.  Sample  Input  files 
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The  following  is  a  sample  Quantum  Espresso  input  file  for  computing  a  converged 
set  of  bands  for  liquid  sodium  using  the  pw.x  executable  from  the  Quantum 
Espresso  suite.  Keyword  descriptions  and  input  syntax  are  available  in  the  Quantum 
Espresso  user  manual. 


&CONTROL 
calculation  =  "scf", 

pseudo_dir  =  "/usr/people/detaylor/PSEUDOPOTENTIALS/ultra", 
prefix  =  "dtpoly", 
nstep  =  1 

wf_collect  =  .true. 

/ 

&SYSTEM 
ibrav  =  0, 
nat  =  54, 
ntyp  =  1, 
ecutwfc  =  30.D0, 
vdw_corr  =  "grimme-d2" 
nbnd=280 

occupations  =  "smearing" 
smearing  =  "fd" 
degauss  =  0.006 

/ 

&ELECTRONS 
electron_maxstep  =  500, 
conv_thr  =  l.D-6, 
mixing_beta  =  0.3D0, 
scf_must_converge  =  .false. 

/ 

&IONS 

/ 

&CELL 

/ 

CELL_PARAMETERS  angstrom 
12.480  0  0 
0  12.480  0 
0  0  12.480 
ATOMIC_SPECIES 
Na  1.00  na.upf 

AT OMIC_POS ITION S  {angstrom} 

Na  -6.49145  -5.56509  7.46734 
Na  -2.63957  9.48177  0.83205 
Na  9.14960  -1.19936  -0.61690 
Na  -2.88286  18.74158  5.66213 
Na  -13.51095  15.30019  13.18334 
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Na 

12.53685 

-3.57070 

9.70781 

Na 

-9.29464 

14.31882 

-8.44320 

Na 

-12.00947 

-4.46376 

0.09270 

Na 

1.70977 

3.09272 

5.64892 

Na 

7.78536 

0.94219 

18.73531 

Na 

-3.50624 

-7.49415 

16.25346 

Na 

2.02219 

-2.22495 

4.69396 

Na 

-6.14664 

7.98090 

4.55669 

Na 

-0.81211 

18.41742 

-15.53245 

Na 

-2.75103 

15.48755 

-0.72738 

Na 

4.25365 

7.57677 

12.72076 

Na 

15.75801 

14.17200 

-1.23555 

Na 

7.68093 

15.46983 

20.90631 

Na 

7.20992 

9.10673 

13.03992 

Na 

18.22595 

0.77448 

-3.65508 

Na 

11.58326 

-0.24269 

10.80184 

Na 

-2.48933  ■ 

-11.93361 

8.67884 

Na 

1.99670- 

-13.11942 

14.66121 

Na 

2.35281 

-1.45250 

11.17339 

Na 

-2.25883 

3.15010 

4.63930 

Na 

4.85383 

5.35284 

-2.36860 

Na 

12.83331 

5.24590 

3.13106 

Na 

7.02326 

4.79217 

-1.14337 

Na 

4.00584 

0.78262 

18.36582 

Na 

12.50268 

12.51815 

4.52492 

Na 

5.78203 

-2.72011 

-5.59039 

Na 

-4.00686 

6.59546 

1.39452 

Na 

12.01031 

11.46302 

1.91613 

Na 

6.74905 

3.14472 

4.42347 

Na 

2.44556 

19.70637 

20.04626 

Na 

17.92583 

15.63488 

7.49345 

Na 

-1.56790 

-2.23829 

7.54851 

Na 

5.69915 

13.74492 

2.27491 

Na 

12.09839 

-3.83016 

5.82794 

Na 

10.72698 

4.07149 

8.18170 

Na 

7.94006 

0.98534 

10.44330 

Na 

19.90926 

-3.95616 

10.43586 

Na 

0.01664 

7.96829 

-9.70577 

Na 

3.30370 

8.37723 

9.82259 

Na 

14.48450 

5.28234 

-1.30440 

Na 

21.55710 

10.20285 

-7.61440 

Na 

13.51033 

4.36419 

-3.36310 

Na 

4.09043 

17.01501 

1.34137 

Na 

2.62843 

7.74397 

4.58241 

Na 

26.78942 

0.63389 

8.15922 

Na 

5.67407 

-1.12595 

-0.10256 
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Na  6.28368  10.45636  15.06059 

Na  4.32227  10.90205  9.00709 

Na  3.73724  4.93973  4.76614 

K_POINTS  {automatic} 

222000 

Below  is  a  sample  Quantum  Espresso  input  file  for  computing  momentum  matrix 
elements  using  the  modified  bands. x  executable.  Keyword  descriptions  and  input 
syntax  are  available  in  the  Quantum  Espresso  user  manual. 


&BANDS 
prefix  =  "dtpoly" 

outdir  =  '7p/workl/detaylor/pw2gw/today/conduct/sodium-liquid" 
filband=  "bands.dat" 
lp  =  .true. 

filp  =  "carlos.mom" 
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List  of  Symbols,  Acronyms,  and  Abbreviations 


ARL  US  Army  Research  Laboratory 

DC  direct  current 

DSI  Director’s  Strategic  Initiative 

eV  electronvolt 

KG  Kubo-Greenwood 

SCF  self-consistent  field 
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